home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / geom_lib / geomat3d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  22.4 KB  |  683 lines
  1. /******************************************************************************
  2. * GeoMat3d.c - Trans. Matrices , Vector computation, and Comp.geom.          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March 1990.                          *
  5. ******************************************************************************/
  6.  
  7. #include <math.h>
  8. #include <stdio.h>
  9. #include "irit_sm.h"
  10. #include "iritprsr.h"
  11. #include "allocate.h"
  12. #include "convex.h"
  13. #include "geomat3d.h"
  14.  
  15. /* #define DEBUG1       Print information on entry and exit of routines. */
  16.  
  17. static int CGPointRayRelation(PointType Pt, PointType PtRay, int Axes);
  18.  
  19. /*****************************************************************************
  20. *  Routine to copy one vector to another:                     *
  21. *****************************************************************************/
  22. void GMVecCopy(VectorType Vdst, VectorType Vsrc)
  23. {
  24.     int i;
  25.  
  26.     for (i = 0; i < 3; i++)
  27.     Vdst[i] = Vsrc[i];
  28. }
  29.  
  30. /*****************************************************************************
  31. *  Routine to normalize the vector length to a unit length:                  *
  32. *****************************************************************************/
  33. void GMVecNormalize(VectorType V)
  34. {
  35.     int i;
  36.     RealType Len;
  37.  
  38.     Len = GMVecLength(V);
  39.     if (Len > 0.0)
  40.         for (i = 0; i < 3; i++) {
  41.         V[i] /= Len;
  42.         if (ABS(V[i]) < EPSILON)
  43.         V[i] = 0.0;
  44.     }
  45. }
  46.  
  47. /*****************************************************************************
  48. *  Routine to calculate the magnitude (length) of a given 3D vector:         *
  49. *****************************************************************************/
  50. RealType GMVecLength(VectorType V)
  51. {
  52.     return sqrt(SQR(V[0]) + SQR(V[1]) + SQR(V[2]));
  53. }
  54.  
  55. /*****************************************************************************
  56. *  Routine to calculate the cross product of two vectors:                    *
  57. * Note Vres might be the same as V1 or V2 !                                  *
  58. *****************************************************************************/
  59. void GMVecCrossProd(VectorType Vres, VectorType V1, VectorType V2)
  60. {
  61.     VectorType Vtemp;
  62.  
  63.     Vtemp[0] = V1[1] * V2[2] - V2[1] * V1[2];
  64.     Vtemp[1] = V1[2] * V2[0] - V2[2] * V1[0];
  65.     Vtemp[2] = V1[0] * V2[1] - V2[0] * V1[1];
  66.  
  67.     GMVecCopy(Vres, Vtemp);
  68. }
  69.  
  70. /*****************************************************************************
  71. *  Routine to calculate the dot product of two vectors:                      *
  72. *****************************************************************************/
  73. RealType GMVecDotProd(VectorType V1, VectorType V2)
  74. {
  75.     return  V1[0] * V2[0] + V1[1] * V2[1] + V1[2] * V2[2];
  76. }
  77.  
  78. /*****************************************************************************
  79. * Routine to generate rotation object around the X axis in Degree degrees:   *
  80. *****************************************************************************/
  81. IPObjectStruct *GMGenMatObjectRotX(RealType *Degree)
  82. {
  83.     MatrixType Mat;
  84.  
  85.     MatGenMatRotX1(DEG2RAD(*Degree), Mat);    /* Generate the trans. matrix. */
  86.  
  87.     return GenMATObject(Mat);
  88. }
  89.  
  90. /*****************************************************************************
  91. * Routine to generate rotation object around the Y axis in Degree degrees:   *
  92. *****************************************************************************/
  93. IPObjectStruct *GMGenMatObjectRotY(RealType *Degree)
  94. {
  95.     MatrixType Mat;
  96.  
  97.     MatGenMatRotY1(DEG2RAD(*Degree), Mat);    /* Generate the trans. matrix. */
  98.  
  99.     return GenMATObject(Mat);
  100. }
  101.  
  102. /*****************************************************************************
  103. * Routine to generate rotation object around the Z axis in Degree degrees:   *
  104. *****************************************************************************/
  105. IPObjectStruct *GMGenMatObjectRotZ(RealType *Degree)
  106. {
  107.     MatrixType Mat;
  108.  
  109.     MatGenMatRotZ1(DEG2RAD(*Degree), Mat);    /* Generate the trans. matrix. */
  110.  
  111.     return GenMATObject(Mat);
  112. }
  113.  
  114. /*****************************************************************************
  115. * Routine to generate translation object according to XYZ vector Vec.         *
  116. *****************************************************************************/
  117. IPObjectStruct *GMGenMatObjectTrans(VectorType Vec)
  118. {
  119.     MatrixType Mat;
  120.  
  121.     /* Generate the transformation matrix */
  122.     MatGenMatTrans(Vec[0], Vec[1], Vec[2], Mat);
  123.  
  124.     return GenMATObject(Mat);
  125. }
  126.  
  127. /*****************************************************************************
  128. * Routine to scale an object according to XYZ scaling vector Vec.         *
  129. *****************************************************************************/
  130. IPObjectStruct *GMGenMatObjectScale(VectorType Vec)
  131. {
  132.     MatrixType Mat;
  133.  
  134.     /* Generate the transformation matrix */
  135.     MatGenMatScale(Vec[0], Vec[1], Vec[2], Mat);
  136.  
  137.     return GenMATObject(Mat);
  138. }
  139.  
  140. /*****************************************************************************
  141. * Routine to transform an object according to the transformation matrix.     *
  142. *****************************************************************************/
  143. IPObjectStruct *GMTransformObject(IPObjectStruct *PObj, MatrixType Mat)
  144. {
  145.     int i;
  146.     IPObjectStruct *NewPObj;
  147.  
  148.     if (IP_IS_POLY_OBJ(PObj)) {
  149.         int IsPolygon = !IP_IS_POLYLINE_OBJ(PObj);
  150.         VectorType Pin, PTemp;
  151.         IPPolygonStruct *Pl;
  152.         IPVertexStruct *V, *VFirst;
  153.  
  154.         NewPObj = CopyObject(NULL, PObj, FALSE);
  155.     Pl = NewPObj -> U.Pl;
  156.  
  157.     while (Pl != NULL) {
  158.         V = VFirst = Pl -> PVertex;
  159.         PT_ADD(Pin, V -> Coord, Pl -> Plane);  /* Prepare point IN side. */
  160.         MatMultVecby4by4(Pin, Pin, Mat); /* Transformed relative to new. */
  161.  
  162.         do {
  163.         if (IsPolygon)
  164.             PT_ADD(PTemp, V -> Coord, V -> Normal);
  165.  
  166.         MatMultVecby4by4(V -> Coord, V -> Coord, Mat);/* Update pos. */
  167.  
  168.         if (IsPolygon) {
  169.             MatMultVecby4by4(PTemp, PTemp, Mat);   /* Update normal. */
  170.             PT_SUB(V -> Normal, PTemp, V -> Coord);
  171.             PT_NORMALIZE(V -> Normal);
  172.         }
  173.  
  174.         V = V -> Pnext;
  175.         }
  176.         while (V != VFirst && V != NULL);           /* VList is circular! */
  177.  
  178.         if (IsPolygon)
  179.         IritPrsrUpdatePolyPlane2(Pl, Pin);    /* Update plane eqn. */
  180.  
  181.         Pl = Pl -> Pnext;
  182.     }
  183.     }
  184.     else if (IP_IS_SRF_OBJ(PObj)) {
  185.     CagdSrfStruct *Srf;
  186.  
  187.         NewPObj = CopyObject(NULL, PObj, FALSE);
  188.  
  189.     for (Srf = NewPObj -> U.Srfs; Srf != NULL; Srf = Srf -> Pnext)
  190.         CagdSrfMatTransform(Srf, Mat);
  191.     }
  192.     else if (IP_IS_CRV_OBJ(PObj)) {
  193.     CagdCrvStruct *Crv;
  194.  
  195.         NewPObj = CopyObject(NULL, PObj, FALSE);
  196.  
  197.     for (Crv = NewPObj -> U.Crvs; Crv != NULL; Crv = Crv -> Pnext)
  198.         CagdCrvMatTransform(Crv, Mat);
  199.     }
  200.     else if (IP_IS_POINT_OBJ(PObj)) {
  201.         NewPObj = CopyObject(NULL, PObj, FALSE);
  202.  
  203.     MatMultVecby4by4(NewPObj -> U.Pt, NewPObj -> U.Pt, Mat);
  204.     }
  205.     else if (IP_IS_VEC_OBJ(PObj)) {
  206.         NewPObj = CopyObject(NULL, PObj, FALSE);
  207.  
  208.     MatMultVecby4by4(NewPObj -> U.Vec, NewPObj -> U.Vec, Mat);
  209.     }
  210.     else if (IP_IS_OLST_OBJ(PObj)) {
  211.     IPObjectStruct *PObjTmp;
  212.  
  213.     NewPObj = IPAllocObject("", IP_OBJ_LIST_OBJ, NULL);
  214.  
  215.         for (i = 0; (PObjTmp = ListObjectGet(PObj, i)) != NULL; i++)
  216.             ListObjectInsert(NewPObj, i, GMTransformObject(PObjTmp, Mat));
  217.     ListObjectInsert(NewPObj, i, NULL);
  218.     }
  219.     else {
  220.         NewPObj = CopyObject(NULL, PObj, FALSE);
  221.     }
  222.  
  223.     return NewPObj;
  224. }
  225.  
  226. /*****************************************************************************
  227. * Routine to transform an object list according to transformation matrix.    *
  228. *****************************************************************************/
  229. IPObjectStruct *GMTransformObjectList(IPObjectStruct *PObj, MatrixType Mat)
  230. {
  231.     IPObjectStruct
  232.     *PTailObj = NULL,
  233.     *PTransObj = NULL;
  234.  
  235.     for ( ; PObj != NULL; PObj = PObj -> Pnext) {
  236.     if (PTailObj == NULL)
  237.         PTailObj = PTransObj = GMTransformObject(PObj, Mat);
  238.     else {
  239.         PTailObj -> Pnext = GMTransformObject(PObj, Mat);
  240.         PTailObj = PTailObj -> Pnext;
  241.     }
  242.     }
  243.  
  244.     return PTransObj;
  245. }
  246.  
  247. /*****************************************************************************
  248. *   Routine to calc the distance between two 3d points                       *
  249. *****************************************************************************/
  250. RealType CGDistPointPoint(PointType P1, PointType P2)
  251. {
  252.     RealType t;
  253.  
  254. #ifdef DEBUG1
  255.     printf("CGDistPointPoint entered\n");
  256. #endif /* DEBUG1 */
  257.  
  258.     t = sqrt(SQR(P2[0]-P1[0]) + SQR(P2[1]-P1[1]) + SQR(P2[2]-P1[2]));
  259.  
  260. #ifdef DEBUG1
  261.     printf("CGDistPointPoint exit\n");
  262. #endif /* DEBUG1 */
  263.  
  264.      return t;
  265. }
  266.  
  267. /*****************************************************************************
  268. *   Routine to create the plane from given 3 points. if two of the points    *
  269. * are the same it returns FALSE, otherwise (succesfull) returns TRUE.         *
  270. *****************************************************************************/
  271. int CGPlaneFrom3Points(PlaneType Plane, PointType Pt1, PointType Pt2,
  272.                             PointType Pt3)
  273. {
  274.     VectorType V1, V2;
  275.  
  276. #ifdef DEBUG1
  277.     printf("CGPlaneFrom3Points entered\n");
  278. #endif /* DEBUG1 */
  279.  
  280.     if (PT_APX_EQ(Pt1, Pt2) || PT_APX_EQ(Pt2, Pt3) || PT_APX_EQ(Pt1, Pt3))
  281.     return FALSE;
  282.  
  283.     PT_SUB(V1, Pt2, Pt1);
  284.     PT_SUB(V2, Pt3, Pt2);
  285.     GMVecCrossProd(Plane, V1, V2);
  286.     PT_NORMALIZE(Plane);
  287.  
  288.     Plane[3] = -DOT_PROD(Plane, Pt1);
  289.  
  290. #ifdef DEBUG1
  291.     printf("CGPlaneFrom3Points exit\n");
  292. #endif /* DEBUG1 */
  293.  
  294.     return TRUE;
  295. }
  296.  
  297. /*****************************************************************************
  298. *   Routine to calc the closest 3d point to a given 3d line. the line is     *
  299. * given as a point on it (Pl) and vector (Vl).                               *
  300. *****************************************************************************/
  301. void CGPointFromPointLine(PointType Point, PointType Pl, PointType Vl,
  302.                             PointType ClosestPoint)
  303. {
  304.     int i;
  305.     PointType V1, V2;
  306.     RealType CosAlfa, VecMag;
  307.  
  308. #ifdef DEBUG1
  309.     printf("CGPointFromLinePlane entered\n");
  310. #endif /* DEBUG1 */
  311.  
  312.     for (i = 0; i < 3; i++) {
  313.         V1[i] = Point[i] - Pl[i];
  314.         V2[i] = Vl[i];
  315.     }
  316.     VecMag = GMVecLength(V1);
  317.     GMVecNormalize(V1);/* Normalized vector from Point to a point on line Pl.*/
  318.     GMVecNormalize(V2);            /* Normalized line direction vector. */
  319.  
  320.     CosAlfa = GMVecDotProd(V1, V2);/* Find the angle between the two vectors.*/
  321.  
  322.     /* Find P1 - the closest point to Point on the line: */
  323.     for (i = 0; i < 3; i++)
  324.     ClosestPoint[i] = Pl[i] + V2[i] * CosAlfa * VecMag;
  325.  
  326. #ifdef DEBUG1
  327.     printf("CGPointFromLinePlane exit\n");
  328. #endif /* DEBUG1 */
  329. }
  330.  
  331. /*****************************************************************************
  332. *   Routine to calc the distance between 3d point and 3d line. the line is   *
  333. * given as a point on it (Pl) and vector (Vl).                               *
  334. *****************************************************************************/
  335. RealType CGDistPointLine(PointType Point, PointType Pl, PointType Vl)
  336. {
  337.     RealType t;
  338.     PointType Ptemp;
  339.  
  340. #ifdef DEBUG1
  341.     printf("CGDistPointLine entered\n");
  342. #endif /* DEBUG1 */
  343.  
  344.     CGPointFromPointLine(Point, Pl, Vl, Ptemp);/* Find closest point on line.*/
  345.     t = CGDistPointPoint(Point, Ptemp);
  346.  
  347. #ifdef DEBUG1
  348.     printf("CGDistPointLine exit\n");
  349. #endif /* DEBUG1 */
  350.  
  351.     return t;
  352. }
  353.  
  354. /*****************************************************************************
  355. *   Routine to calc the distance between a Point and a Plane. The Plane is   *
  356. * defined with 4 coef. : Ax + By + Cz + D = 0 given as 4 elements vector.    *
  357. *****************************************************************************/
  358. RealType CGDistPointPlane(PointType Point, PlaneType Plane)
  359. {
  360.     RealType t;
  361.  
  362. #ifdef DEBUG1
  363.     printf("CGDistPointPlane entered\n");
  364. #endif /* DEBUG1 */
  365.  
  366.     t = ((Plane[0] * Point [0] +
  367.       Plane[1] * Point [1] +
  368.       Plane[2] * Point [2] +
  369.       Plane[3]) /
  370.      sqrt(SQR(Plane[0]) + SQR(Plane[1]) + SQR(Plane[2])));
  371.  
  372. #ifdef DEBUG1
  373.     printf("CGDistPointPlane exit\n");
  374. #endif /* DEBUG1 */
  375.  
  376.     return t;
  377. }
  378.  
  379. /*****************************************************************************
  380. *   Routine to find the intersection point of a line and a plane (if any)    *
  381. *   The Plane is defined with 4 coef. : Ax + By + Cz + D = 0 given as 4      *
  382. * elements vector. The line is define via a point on it Pl and a direction   *
  383. * vector Vl. Return TRUE only if such point exists.                          *
  384. *****************************************************************************/
  385. int CGPointFromLinePlane(PointType Pl, PointType Vl, PlaneType Plane,
  386.                     PointType InterPoint, RealType *t)
  387. {
  388.     int i;
  389.     RealType DotProd;
  390.  
  391. #ifdef DEBUG1
  392.     printf("CGPointFromLinePlane entered\n");
  393. #endif /* DEBUG1 */
  394.  
  395.     /* Check to see if they are vertical - no intersection at all! */
  396.     DotProd = GMVecDotProd(Vl, Plane);
  397.     if (ABS(DotProd) < EPSILON)
  398.     return FALSE;
  399.  
  400.     /* Else find t in line such that the plane equation plane is satisfied: */
  401.     *t = (-Plane[3] - Plane[0] * Pl[0] - Plane[1] * Pl[1] - Plane[2] * Pl[2])
  402.                                 / DotProd;
  403.  
  404.     /* And so find the intersection point which is at that t: */
  405.     for (i = 0; i < 3; i++)
  406.     InterPoint[i] = Pl[i] + *t * Vl[i];
  407.  
  408. #ifdef DEBUG1
  409.     printf("CGPointFromLinePlane exit\n");
  410. #endif /* DEBUG1 */
  411.  
  412.     return TRUE;
  413. }
  414.  
  415. /*****************************************************************************
  416. *   Routine to find the intersection point of a line and a plane (if any)    *
  417. *   The Plane is defined with 4 coef. : Ax + By + Cz + D = 0 given as 4      *
  418. * elements vector. The line is define via a point on it Pl and a direction   *
  419. * vector Vl: Line == Pl + Vl * t, where t is the line parameter.         *
  420. *   Return TRUE only if such point exists, in the t parameter range [0..1].  *
  421. *****************************************************************************/
  422. int CGPointFromLinePlane01(PointType Pl, PointType Vl, PlaneType Plane,
  423.                     PointType InterPoint, RealType *t)
  424. {
  425.     int i;
  426.     RealType DotProd;
  427.  
  428. #ifdef DEBUG1
  429.     printf("CGPointFromLinePlane01 entered\n");
  430. #endif /* DEBUG1 */
  431.  
  432.     /* Check to see if they are vertical - no intersection at all! */
  433.     DotProd = GMVecDotProd(Vl, Plane);
  434.     if (ABS(DotProd) < EPSILON)
  435.     return FALSE;
  436.  
  437.     /* Else find t in line such that the plane equation plane is satisfied: */
  438.     *t = (-Plane[3] - Plane[0] * Pl[0] - Plane[1] * Pl[1] - Plane[2] * Pl[2])
  439.                                 / DotProd;
  440.  
  441.     if ((*t < 0.0 && !APX_EQ(*t, 0.0)) ||      /* Not in parameter range. */
  442.     (*t > 1.0 && !APX_EQ(*t, 1.0)))
  443.     return FALSE;
  444.  
  445.     /* And so find the intersection point which is at that t : */
  446.     for (i = 0; i < 3; i++)
  447.     InterPoint[i] = Pl[i] + *t * Vl[i];
  448.  
  449. #ifdef DEBUG1
  450.     printf("CGPointFromLinePlane01 exit\n");
  451. #endif /* DEBUG1 */
  452.  
  453.     return TRUE;
  454. }
  455.  
  456. /*****************************************************************************
  457. *   Routine to find the two point Pti on the lines (Pli, Vli) ,   i = 1, 2   *
  458. * with the minimal euclidian distance between them. In other words the       *
  459. * distance between Pt1 and Pt2 is defined as distance between the two lines. *
  460. *   The two points are calculated using the fact that if V = (Vl1 cross Vl2) *
  461. * then these two points are the intersection point between the following:    *
  462. * point 1 - a plane (defined by V and line1) and the line line2.             *
  463. * point 2 - a plane (defined by V and line2) and the line line1.             *
  464. *   This function returns TRUE iff the two lines are not parallel!           *
  465. *****************************************************************************/
  466. int CG2PointsFromLineLine(PointType Pl1, PointType Vl1,
  467.               PointType Pl2, PointType Vl2,
  468.               PointType Pt1, RealType *t1,
  469.               PointType Pt2, RealType *t2)
  470. {
  471.     int i;
  472.     PointType Vtemp;
  473.     PlaneType Plane1, Plane2;
  474.  
  475. #ifdef DEBUG1
  476.     printf("CG2PointsFromLineLine entered\n");
  477. #endif /* DEBUG1 */
  478.  
  479.     GMVecCrossProd(Vtemp, Vl1, Vl2);   /* Check to see if they are parallel. */
  480.     if (GMVecLength(Vtemp) < EPSILON) {
  481.     for (i = 0; i < 3; i++)
  482.         Pt1[i] = Pl1[i];             /* Pick point on line1 and find */
  483.     CGPointFromPointLine(Pl1, Pl2, Vl2, Pt2); /* closest point on line2. */
  484.         return FALSE;
  485.     }
  486.  
  487.     /* Define the two planes: 1) Vl1, Pl1, Vtemp and 2) Vl2, Pl2, Vtemp         */
  488.     /* Note this sets the first 3 elements A, B, C out of the 4...         */
  489.     GMVecCrossProd(Plane1, Vl1, Vtemp);            /* Find the A, B, C coef.'s. */
  490.     GMVecNormalize(Plane1);
  491.     GMVecCrossProd(Plane2, Vl2, Vtemp);            /* Find the A, B, C coef.'s. */
  492.     GMVecNormalize(Plane2);
  493.  
  494.     /* and now use a point on the plane to find the 4th coef. D: */
  495.     Plane1[3] = (-GMVecDotProd(Plane1, Pl1)); /* VecDotProd uses only first  */
  496.     Plane2[3] = (-GMVecDotProd(Plane2, Pl2)); /* three elements in vec.      */
  497.  
  498.     /* Thats it! now we should solve for the intersection point between a    */
  499.     /* line and a plane but we already familiar with this problem...         */
  500.     i = CGPointFromLinePlane(Pl1, Vl1, Plane2, Pt1, t1) &&
  501.     CGPointFromLinePlane(Pl2, Vl2, Plane1, Pt2, t2);
  502.  
  503. #ifdef DEBUG1
  504.     printf("CG2PointsFromLineLine exit\n");
  505. #endif /* DEBUG1 */
  506.  
  507.     return i;
  508. }
  509.  
  510. /*****************************************************************************
  511. *   Routine to find the distance between two lines (Pli, Vli) ,  i = 1, 2.   *
  512. *****************************************************************************/
  513. RealType CGDistLineLine(PointType Pl1, PointType Vl1,
  514.             PointType Pl2, PointType Vl2)
  515. {
  516.     RealType t1, t2;
  517.     PointType Ptemp1, Ptemp2;
  518.  
  519. #ifdef DEBUG1
  520.     printf("CGDistLineLine entered\n");
  521. #endif /* DEBUG1 */
  522.  
  523.     CG2PointsFromLineLine(Pl1, Vl1, Pl2, Vl2, Ptemp1, &t1, Ptemp2, &t2);
  524.     t1 = CGDistPointPoint(Ptemp1, Ptemp2);
  525.  
  526. #ifdef DEBUG1
  527.     printf("CGDistLineLine exit\n");
  528. #endif /* DEBUG1 */
  529.  
  530.     return t1;
  531. }
  532.  
  533. /*****************************************************************************
  534. *   Routine implements Jordan Theorem: Fire a ray from a given point and find*
  535. * number of intersections of ray with the polygon, excluding the given point *
  536. * Pt (start of ray) itself, if on polygon boundary. The ray is fired in +X   *
  537. * (Axes == 0) or +Y (Axes == 1). And only the X/Y coordinates of the polygon *
  538. * are taken into account, i.e. the orthogonal projection of the polygon on   *
  539. * a X/Y parallel plane (equal to polygon itself if on X/Y parallel plane...) *
  540. *   Note that if the point is on polygon boundary, the ray should not be in  *
  541. * its edge direction                                 *
  542. *                                         *
  543. *   Algorithm:                                     *
  544. * 1. Set NumOfIntersection = 0;                             *
  545. *    Find vertex V not on Ray level and set AlgState to its level (below or  *
  546. *    above the ray level). If none goto 3                     *
  547. *    Mark VStart = V;                                 *
  548. * 2. 2.1. While State(V) == AlgState do                         *
  549. *        2.1.1. V = V -> Pnext;                         *
  550. *        2.1.2. If V == VStart goto 3                     *
  551. *      end;                                     *
  552. *      IntersectionMinX = INFINITY;                         *
  553. *    2.2. While State(V) == ON_RAY do                         *
  554. *        IntersectionMin = MIN(IntersectionMin, V -> Coord[Axes]);    *
  555. *        V = V -> Pnext;                             *
  556. *         end;                                     *
  557. *    2.3. If State(V) != AlgState do                         *
  558. *        2.3.1. Find the intersection point between polygon edge      *
  559. *               Vlast, V and the Ray and update IntersectionMin if    *
  560. *               lower than it.                         *
  561. *        2.3.2. If IntersectionMin is greater than Pt[Axes] increase  *
  562. *               the NumOfIntersection counter by 1.             *
  563. *      end;                                     *
  564. *    2.4. AlgState = State(V);                             *
  565. *    2.5. goto 2.1.                                 *
  566. * 3. return NumOfIntersection;                             *
  567. *                                         *
  568. *****************************************************************************/
  569. int CGPolygonRayInter(IPPolygonStruct *Pl, PointType PtRay, int RayAxes)
  570. {
  571.     int NewState, AlgState, RayOtherAxes,
  572.     NumOfInter = 0,
  573.     Quit = FALSE;
  574.     RealType InterMin, Inter, t;
  575.     IPVertexStruct *V, *VStart,
  576.     *VLast = NULL;
  577.  
  578.     RayOtherAxes = (RayAxes == 1 ? 0 : 1);     /* Other dir: X -> Y, Y -> X. */
  579.  
  580.     /* Stage 1 - find a vertex below the ray level: */
  581.     V = VStart = Pl -> PVertex;
  582.     do {
  583.     if ((AlgState = CGPointRayRelation(V -> Coord, PtRay, RayOtherAxes))
  584.                             != ON_RAY)
  585.         break;
  586.     V = V -> Pnext;
  587.     }
  588.     while (V != VStart && V != NULL);
  589.     if (AlgState == ON_RAY)
  590.     return 0;
  591.     VStart = V; /* Vertex Below Ray level */
  592.  
  593.     /* Stage 2 - scan the vertices and count number of intersections. */
  594.     while (!Quit) {
  595.     /* Stage 2.1. : */
  596.     while (CGPointRayRelation(V -> Coord, PtRay, RayOtherAxes) == AlgState) {
  597.         VLast = V;
  598.         V = V -> Pnext;
  599.         if (V == VStart) {
  600.         Quit = TRUE;
  601.         break;
  602.         }
  603.     }
  604.     InterMin = INFINITY;
  605.  
  606.     /* Stage 2.2. : */
  607.     while (CGPointRayRelation(V -> Coord, PtRay, RayOtherAxes) == ON_RAY) {
  608.         InterMin = MIN(InterMin, V -> Coord[RayAxes]);
  609.         VLast = V;
  610.         V = V -> Pnext;
  611.         if (V == VStart)
  612.         Quit = TRUE;
  613.     }
  614.  
  615.     /* Stage 2.3. : */
  616.     if ((NewState = CGPointRayRelation(V -> Coord, PtRay, RayOtherAxes))
  617.                                 != AlgState) {
  618.         /* Stage 2.3.1 Intersection with ray is in middle of edge: */
  619.         t = (PtRay[RayOtherAxes] - V -> Coord[RayOtherAxes]) /
  620.         (VLast -> Coord[RayOtherAxes] - V -> Coord[RayOtherAxes]);
  621.         Inter = VLast -> Coord[RayAxes] * t +
  622.             V -> Coord[RayAxes] * (1.0 - t);
  623.         InterMin = MIN(InterMin, Inter);
  624.  
  625.         /* Stage 2.3.2. comp. with ray base and inc. # of inter if above.*/
  626.         if (InterMin > PtRay[RayAxes] &&
  627.         !APX_EQ(InterMin, PtRay[RayAxes]))
  628.         NumOfInter++;
  629.     }
  630.  
  631.     AlgState = NewState;
  632.     }
  633.  
  634.     return NumOfInter;
  635. }
  636.  
  637. /*****************************************************************************
  638. *   Routine to return the relation between the ray level and a given point,  *
  639. * to be used in the CGPolygonRayInter routine above.                 *
  640. *****************************************************************************/
  641. static int CGPointRayRelation(PointType Pt, PointType PtRay, int Axes)
  642. {
  643.     if (APX_EQ(PtRay[Axes], Pt[Axes]))
  644.         return ON_RAY;
  645.     else if (PtRay[Axes] < Pt[Axes])
  646.         return ABOVE_RAY;
  647.     else
  648.     return BELOW_RAY;
  649. }
  650.  
  651.  
  652. /*****************************************************************************
  653. * Same as CGPolygonRayInter but for arbitrary oriented polygon.             *
  654. * The polygon (and the point) is first rotated to a XY parallel plane.         *
  655. *****************************************************************************/
  656. int CGPolygonRayInter3D(IPPolygonStruct *Pl, PointType PtRay, int RayAxes)
  657. {
  658.     int i;
  659.     MatrixType RotMat;
  660.     IPVertexStruct *V, *VHead;
  661.     IPPolygonStruct *RotPl;
  662.     PointType RotPt;
  663.  
  664.     /* Make a copy of original to work on. */
  665.     RotPl = IPAllocPolygon(1, Pl ->Tags, CopyVertexList(Pl -> PVertex), NULL);
  666.  
  667.     /* Bring the polygon to a XY parallel plane by rotation. */
  668.     GenRotateMatrix(RotMat, Pl -> Plane);
  669.     V = VHead = RotPl -> PVertex;
  670.     do {                    /* Transform the polygon itself. */
  671.     MatMultVecby4by4(V -> Coord, V -> Coord, RotMat);
  672.     V = V -> Pnext;
  673.     }
  674.     while (V != NULL && V != VHead);
  675.     MatMultVecby4by4(RotPt, PtRay, RotMat);
  676.  
  677.     i = CGPolygonRayInter(RotPl, RotPt, RayAxes);
  678.  
  679.     IPFreePolygonList(RotPl);
  680.  
  681.     return i;
  682. }
  683.